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Abstract 

We introduce a new strategy for coupling the parallel in time (parareal) iterative 
methodology with multiscale integrators. Following the parareal framework, the algorithm 
computes a low-cost approximation of all slow variables in the system using an appropriate 
multiscale integrator, which is refined using parallel fine scale integrations. Convergence is 
obtained using an alignment algorithm for fast phase-like variables. The method may be 
used either to enhance the accuracy and range of applicability of the multiscale method in 
approximating only the slow variables, or to resolve all the state variables. The numerical 
scheme does not require that the system is split into slow and fast coordinates. Moreover, 
the dynamics may involve hidden slow variables, for example, due to resonances. We 
propose an alignment algorithm for almost-periodic solution, in which case convergence 
of the parareal iterations is proved. The applicability of the method is demonstrated in 
numerical examples. 


1 Introduction 

The parallel in time, also known as the “parareal” method, introduced by Lions, Maday and 
Turinici }49| is a simple yet effective scheme for the parallelization of numerical solutions for 
a large class of time dependent problems [50] • It consists of a fixed point iteration involving 
a coarse-but-cheap and a hne-but-expensive integrators. Computational time is reduced by 
parallelization of the fine integrations. For problems with separated multiple scales, it is 
tempting to apply a multiscale solver as a coarse integrator. So far, such types of parallel 
methods are limited to a few special multiscale cases such as chemical kinetics HU23I33], 
dissipative ordinary differential equations (ODEs) m and highly oscillatory (HiOsc) problems 
in which the oscillatory behavior is relatively simple [20 L [32j . One difficulty stems out from 
a fundamental difference between the parareal and the multiscale philosophies — while the 
former requires point-wise convergence of the numerical solvers (in the state variable), most 
multiscale schemes gain efficiency by only approximating a reduced set of slowly varying 
coarse/slow/macroscopic variables [101 [TO '221 [29l fill [55] . 
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In this paper, we develop a general strategy that couples multiscale integrators and fully 
resolved fine scale integration for parallel in time computation of HiOsc solutions of a class of 
ODEs. There are several advantages in such coupling strategies. First, some multiscale meth¬ 
ods (such as the Poincare-map technique [2]) only approximate the slow constituents or slow 
variables of the dynamics. Proper coupling of multiscale and fine scale solvers via a parareal- 
like framework can be efficient (by parallelization) in computing full detailed solutions, in¬ 
cluding the fast phase in the HiOsc dynamics. The choice of multiscale method is not limited 
to the Poincare-map technique. Any multiscale method can be used as a coarse integrator 
as long as fast phase-like variables are appropriately aligned as required. Then, convergence 
of the parareal iterations can be shown in a similar manner. Second, the parareal iterations 
enhance the stability and accuracy of the multiscale scheme, in particular when the scale sep¬ 
aration in the system is not significant and the corresponding sampling/averaging errors are 
non-negligible. Finally, parareal multiscale coupling schemes can deal with more challenging 
situations, for example, (a) the effective equation is valid almost everywhere macroscopically, 
but is not an adequate description of the system at small but a priori ’’unpredictable” locations 
in the phase space (as these regions may depend on the solutions)] and (b) the influence of 
microscopic solutions in these regions on the macroscopic solution elsewhere is significant. 

In @7], Legoll et at suggest a multiscale parareal scheme for singularly perturbed ODEs in 
which the fast dynamics is dissipative, i.e., the dynamics relaxes rapidly to a low dimensional 
manifold. One of the main contributions of m is the understanding that the slow and fast 
parts of the dynamics need to be addressed separately. They suggest two approaches: The first 
is a straight-forward application of parareal, which is shown to converge but loses accuracy as 
the system becomes more singular. In Section fl ,31 we demonstrate that naive parareal does 
not converge when applied for HiOsc systems. The second approach assumes that the system 
is split into slow and fast variables, or that a change of variables splitting the system is given. 
This approach may be applied to HiOsc systems, but it is relatively restrictive as in many 
examples and applications such a splitting is not known. Dai et al m suggest an application 
of the parareal framework to Hamiltonian systems. They consider two main approaches. The 
first is a time-reversible iteration scheme (applied together with time-reversible fine and coarse 
integrators). The second projects solutions at coarse time segments onto the constant energy 
manifold. The two approaches are also combined together. The first approach is specific to 
Hamiltonian dynamics and not to general HiOsc problems. The projection method cannot be 
applied to the HiOsc case because the main difficulty is not with the approximation of slow 
variables (or constants of motion), but with the fast phase. In addition, since the methods 
presented in [20] are not multiscale, their accuracy and efficiency are expected to deteriorate 
when the frequencies of oscillations are large. Combining the symmetric approach of Dai et al 
with our alignment method for Hamiltonian systems may be an interesting application, but is 
beyond the scope of the current manuscript. In particular, the alignment algorithms should 
also be made symmetric, similar to the ideas of Dai et al. Applications of parareal methods 
to Hamiltonian dynamics is also analyzed in [27j. Additional approaches to use symplectic 
integrators with applications to molecular dynamics include pans E8]. Finally, Haut and 
Wingate [32] suggest a parareal method for PDEs with linear HiOsc forcing, As in [37], their 
method applies exact knowledge of the fast variable (the phase in the HiOsc case) to design 
a convergent parareal scheme. In this respect, the method proposed here goes further and is 
also applicable to nonlinear HiOsc systems. However, in this paper the discussion is restricted 
to the ODE case. One of the main goals of the current paper is to design a convergent parareal 
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algorithm that does not require explicit knowledge of the fast and slow variables. 

We begin with a short overview of the parareal method within the context of ODEs and 
test its performance on a simple example HiOsc system. 


1.1 The parareal method for ODEs 

Consider the following initial value problem 

u = f(t,u ), u(0)=u o , (1.1) 

where u € M. d and t £ [0, T], We assume that / is sufficiently smooth. Let H denote an 
intermediate time step, 0 < H < T and N = T/H an integer. Suppose that we are given two 
approximate integrators for (11.11) : a cheap coarse integrator with low accuracy denoted C, and 
a fine, high accuracy integrator which is relatively expensive in terms of efficiency, denoted 
T. The approximate propagation operators to time H obtained using the the coarse and fine 
integrators are denoted by Ch and J~h, respectively. 

Furthermore, denote by u the approximation for u(nH) at the fc’th iteration. For all 
iterations, the initial values are the same VT^Uq = u$. The objective is to have —» F n HU o 
as k — > oo, i.e., convergence to the approximation given by the high-accuracy fine integrator. 
The parareal approximation to (11.11) is as follows. See Figure [T] for a sketch of the parareal 
methodology. 

Algorithm 1.1. 


1. Initialization: Construct the zero’th iteration approximation using a chosen coarse in¬ 
tegrator: 


Uq = uq and = Chu^_i, n = 1,..., N. 


2. Iterations: k = 1... K 

Uq = u 0 and u k n = C H u^-i + F H u k n Z\ - C H u^z\, n = l,...,N. (1.2) 

Note that the calculation of the fine integrator Fhu^z\ in m requires only the initial 
condition which depends on the previous iteration. Hence, for each k, Ftvf^zX, 0 < t < 

H, n = 1, 2, • • ■ ,N can be computed in parallel. The solution computed by the accurate but 
expensive integrator is a fixed point. Indeed, when the iteration is sufficiently large (k >n), 
the solution u^ become identical to it: 

Un = F nH uo, n < k. 

In fact, (11.21) can be regarded as a fixed-point iteration. In [50j, it is proved that under some 
sufficient conditions of /, which we shall recall in Section fl.2l 

\u*-u(nH)\<C(H k + E f ), (1.3) 

where Ef is the global error in solving the full ODE using the fine propagator, and C depends 
on the derivatives of the solutions. Eq. (11.31) assumes a 1st order coarse integrator. 

In order to identify the source of the difficulty in developing parareal algorithms for highly 
oscillatory problems, we adapt the parareal proof of convergence given by Maday in [50] for 
non-singular ODEs. 
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C > 0 


Figure 1: A sketch depicting the par areal methodology. Each parareal iteration is constructed 
using three integrations: Fine integration starting at u^ -1 , coarse integration starting at 
and another coarse integration starting at u^. The first two depend only on previous iterations 
and can therefore be computed in parallel. 

1.2 Convergence of parareal 

We consider ODEs of the form <o> with initial conditions w(0) = uq € D C M rf . We are 
interested in solving (11.11) in a fixed time segment [0, T], The solution is denoted u(t\ uq), t G 
[0, T], Let <F denote the flow map (propagator) associated with (II.ID . 

= u(t\x ), Vt > 0. 

For sufficiently smooth / we have that |<3?fX — <&ty\ < e ct \x — y |. In the following C denotes 
a generic positive constant which may depend on T. Since t < T, the prefactor e ct can be 
bounded by 1 + T~ 1 e CT t. This yields a linear stability bound for !>t, 

\$tx - < (1 + Ct)\x - y\. 

For simplicity, we assume that the coarse integrator Ct is a one-step method with step 
size H while the fine integrator J~t has step size h <C H. In addition, we make the following 
accuracy and stability assumptions on the numerical integrators: 

\Ttx — § t x\ < CtEf(l + |x|), \C t x — ® t x\ < CtE c (l + \x\) (1.4) 

where Ef and E c denote the global sup error in solving m in [0,T] using respectively the 
fine and coarse integrators in the entire domain of interest D. Note that both Ef and E c 
typically depend on T. In addition, 

\E t x — E t y\ < (1 + tC)\x — y\, \C t x - C t y\ < (1 + tC)\x - y\ (1.5) 

Let 5Ft = <h t — Ft and 5Ct = 4> f — Ct, denote the errors in the fine and coarse propagators, 
respectively. Then, by a triangle inequality, 

\5F t x - 5F t y\ < (1+ tCE f )\x - y\, \5C t x - SC t y\ < (1 + tCE c )\x - y\. (1.6) 
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We recite the following theorem from [50]. 

Theorem 1.2. Let K < N/2 = T/2H. Then, for all k < K, 

sup | u k - E nH u 0 \ < C(E c ) k . 

n=0,...,N 


Consequently, 


sup | u k n - u(nH) | < C ( E c ) k + E f 

n=0,...,N 1 


(1.7) 


Proof: Applying the parareal iterations (11.21) . 


U n — EnHUo = 


ChU^-I ~ CHUn-1 + ^HU k n _\ - Eh^u-^HUQ 


ChW^-I ~ - &Ch {fE{n-l)H^o) ~ & C H U n _ L 


( 1 . 8 ) 


SEnu n _\ - 5 Eh (^"{n-i )h u o) 


Using assumption (11.61) . and denoting 9 k = C( 1 + CH) k n {Ej + E c ) k H k | u!f — E n nUo\ , we 
have 9 k < 9 k _ x + 9 k z\ ■ By induction, 9 k < C ^ ^ ■ Assuming that Ef < E c , 


u n ~ FnHUQ 


< CT K E k = 0(E; 


(1.9) 


1.3 Parareal and HiOsc ODEs 

We consider HiOsc ODEs given in the singular perturbation form 

u = e~ 1 f 1 (u) +f 0 (u), ( 1 . 10 ) 

with initial condition u(0) = uq € D C M rf , where D is a domain uniformly bounded in e. 
The parameter 0 < e < eo < 1 characterizes the separation of time scales - the fast scale 
involves oscillations with frequencies of order e -1 while the computational time domain is 
[0, T] with T independent of e. Throughout the paper we assume that /i, /o are sufficiently 
smooth, and that for each uq € D, u(t ) is uniformly bounded in e in the time interval [0, T]. 
Furthermore, we assume that the Jacobian of f\ has only purely imaginary eigenvalues in D, 
which are bounded away from 0 and independent of e. These settings typically imply that 
the computational complexity of direct non-multiscale methods is at least ©(e -1 ). 

To understand some of the challenges in applying the parareal framework to HiOsc sys¬ 
tems, we consider the following simple example 

u = (a + ie~ 1 )u, u(0) = 1. (1.11) 

With a > 0, the trajectory of u(t) = e^ a+te ^ is an expanding spiral in the complex plane. 
We further assume that the fine integrator is exact, EfU = e^ a+ie ^u. We first investigate the 
performance of Algorithm II.II using two conventional methods as Ct- Implicit Euler, Explicit 
Euler, and Trapezoidal Rule. Table Q] compares the minimal number of parareal iterations, 
K , to reach an absolute error below 1/10. We observe that when conventional methods are 
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Coarse integrator 

e= 

0.2 

0.1 

0.05 

0.02 

0.01 

0.001 

Explicit Euler (H = e/5) 

K 

7 

12 

22 

52 

607 

12200 

Explicit Euler (H = 1/10) 

K 

34 

79 

100 

100 

100 

100 

Implicit Euler (H = e/5) 

K 

6 

8 

13 

25 

44 

351 

Implicit Euler (H = 1/10) 

K 

18 

49 

93 

100 

100 

100 

Trapezoidal Rule (H = e/5) 

K 

1 

1 

2 

3 

5 

29 

Trapezoidal Rule (H = 1/10) 

K 

4 

18 

71 

100 

100 

100 

The proposed method (H = 1/10) 

K 

1 

1 

1 

1 

1 

1 


Table 1: The number of parareal iterations required to yield an absolute errors of 1/10 in 
the expanding spiral example. Parameters are a = 1/10, T = 10. The maximal number of 
iterations is T/H. 


implemented as a coarse integrator K becomes prohibitively large as e gets small. The increase 
in K for conventional coarse integrators can be explained by the error estimate (11.311 . The 
difficulty lies in the constant C, which grows rapidly with 1/e. For an order p coarse integrator, 
the error is proportional to the p + 1 time derivative of /, which is of order 0(e~( p+1 )). As a 
result, the parareal error for HiOsc systems (HID depends on e, 

I < - u(nH) | <c{E f + [e- 1 (e" 1 ^]*} . (1.12) 

An immediate consequence is that H has to be o(e), even when applying A-stable or symplectic 
methods. See for example the conclusion in (20) . 

This simple example reveals the reason why a naive implementation of the parareal ap¬ 
proach may not be effective for integrating HiOsc problems: both stability and accuracy 
restrictions require that the coarse integrator take steps of order e. As a result, the number of 
coarse steps is 0(e^ 1 ) and the method may take 0(e _1 ) iterations to converge. For compar¬ 
ison, we also include in Table Q] the results obtained using the proposed multiscale parareal 
method. 

1.4 Layout 

The layout of the paper is as follows. Given a conventional parareal method in (|1.2I) . 

u k n = Chu^-! + Thv^X - C H u k z\, 

Section [2] presents the main difficulty in using a multiscale method as the coarse integrator in 
the parareal framework. Section [3] suggests a general approach for overcoming this difficulty. 
In particular, two versions of the update using the fine solution will be proposed: 

Algorithm 13.21 Jacobi style update for approximating slow variables, 

u n = C [ [U k _ 1 4- F H u^Z\ - C H u k z\ 

v\ = \ 

u k = S 0 {M H u^ l _ 1 ]v i) +ui - So(M H u k Z\]v i) 
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Algorithm 13.51 Gauss-Seidel style update for approximating a state variable, 


u n = C hu ^_i + T h u*_\ - C H u*_ 


,k -1 


vi=S^{F H utz\-^z\,u*) 


u* = So(Mhu*_ i; ui) + Vi - S 0 {M h vi]v 1 ) 



where u* is an approximation of u((n — 1 )H) from the previous step. Section 0] describes the 
implementations of this strategy: 5o(uo; vo) - local alignment in Section (4~T1 and S^(ui]Uq,vq) 
- forward alignment in Section 14.21 Section [5] reviews the Poincare method, a multiscale 
numerical method for efficient integration of HiOsc ODEs presented in [2]. This method will 
be used as a coarse solver in the numerical examples presented in Section [6] We conclude in 
Section [7] 


2 Fast oscillations and parareal 

In order to facilitate the presentation of the main algorithms, we shall first describe the setting 
for the underlying multiscale methods. 

The literature on efficient numerical integration of problems with separated time scale is 
rapidly growing. For HiOsc ODEs, recent approaches include envelope methods [52], FLow 
Averaging integratORS [55] , Young measure mm and equation free approaches m, Mag¬ 
nus methods DUE], Filon methods [36];32]) spectral methods [35y.48], asymptotic expansions 
[HEZj and the Heterogeneous Multiscale Methods m eh [ 22 ]. For a recent review see [23]. 

Typically, multiscale methods tackle the computational difficulty in solving HiOsc ODEs 
by taking advantage of scale separation, and aim at computing only the slowly varying prop¬ 
erties of the oscillatory solutions. It requires that enough information about the influence of 
fast scales on the slower scale dynamics can be obtained by performing localized simulations 
over short times, and thereby better efficiency is achieved. The numerical complexity of these 
methods is therefore much smaller than direct simulations of the given systems with HiOsc 
solutions. For example, [5] presents multiscale algorithms that compute the effective behavior 
of HiOsc dynamical systems by using slow variables that are predetermined either analytically 
or numerically. More precisely, we define a slow variable for the system (11.1011 with solution 
u{t] e) as follows. 

Definition 2.1. A smooth function a(t,e) is called slow to order v > 1 if \d v a/dt u \ < C in 
t € [0,T] for some constants C and T independent of e € (0,eo], eo > 0. A smooth function 
£(u) : D — >• M is called a slow variable with respect to u(t ) if f(t) = £(u(f;e)) is slow to 
order 1. 

See [6] m EH E3 [28| 04] 05] 06] for similar definitions and applications for HiOsc problems. 

In this paper, we will work with the following main assumption. 
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Assumption 2.2. There exists a diffeomorphism '1 1 : u —>■ (£(w), <j>(u)), independent of e, 
separating slow and fast variables such that (£, (j>) along the trajectories of (11.101) satisfies an 
ODE of the form 


£ = 0o(£,</>), £(0) = £(«o), 

0 = e _1 5i(C) +92 (£,<£)> ^*(0) = 4>( u o), 


where f € M d_n , <f> € M n , and 0 < e < £o < 1 is a small parameter. We assume that for 
fixed slow coordinates the fast variable f is ergodic 0 with respect to an invariant manifold, 
which is diffeomorphic to a n-dimensional torus, T n . 


Using ergodicity, one can invoke a theory of averaging [53], which implies that the dynamics 
of slow variables can be approximated (0(e) in the sup norm for 0 < t < T, T = 0(1)) by an 
averaged equation of the form 


I = F(fi), F{£) = J g 0 (£,«£)<% 

1(0) = £(«o)i 


where denotes the invariant measure for cf> at fixed For example, perturbed integrable 
Hamiltonian systems constitute a wide class of dynamical systems that satisfy this assumption. 
From now on, we shall refer to (j> as the phase of u. 

The main objective of many multiscale methods is efficient numerical approximations of 
£(u(t)) only. The general strategy of our algorithm is based on such multiscale methods for 
HiOsc ODEs that only resolve the macroscopic behavior of a system as specified by the slow 
variables [H 0) ;5] IS [7, 8[ [25, o4] . In this respect, the algorithms listed above are different 
from other multiscale methods that resolve all scales of the dynamics, for example, multi-level 
methods or high-order asymptotic expansions [m ns ns sa sa¬ 
lt is possible to design a parareal algorithm for computing only the averaged slow variables 
using multiscale integrators as both the coarse and fine integrators. Such an approach is 
essentially a parareal scheme for the averaged equation. However, this is not the point of 
this paper — here we are interested in the possibility of creating a parareal algorithm that 
computes all state variables, including the fast phase information. 

We consider the problem of using a multiscale integrator in the coarse integration, and 
provide the stability of the corresponding coupling of nrultiscale-fine integrators under the 
parareal framework. Since the error bound stated in (|1.12l) still formally applies in this case, 
one cannot expect convergence of u(t) unless some additional improvement is made to the 
chosen existing multiscale scheme. 

Consider again the simple expanding spiral (11.111) with a = 1. It is easily verified 
that \u(t)\ = e* is a slow variable. For convenience of the discussion, we assume that the 
fine/microscopic solver is exact, i.e. J~tu = e^ 1+ie ^u, and that the coarse/macroscopic 
solver is exact in the slow variables, i.e. any function of \u\ is computed without error but 
the phase of u may be wrong. We write the macroscopic solution as CfiJ = e t e^ tt +9t ' ) U, 
where 0t € [0, 27r) denotes the error in the phase that is produced by the macroscopic solver. 


x By ergodic, we mean that any trajectory of 0(f) can get arbitrarily close to any point in the invariant 
manifold. In particular, this implies the existence of a unique invariant distribution and Birkhoff’s ergodic 
theorem. 









Applying Algorithm 11.11 we obtain 

4 1} = «(3tf) (1 + O(0 2 h )) , 4 2) = «(3fT)(l + C7(^)). (2.3) 

This simple exercise shows that the naive iterations improve the accuracy of the macroscopic 
solution if Oh is small, and that the iterations diverge if Oh is not sufficiently small. However, 
in a typical HiOsc, Oh is not necessarily small. In general, Oh can be any value in [0, 2n\ and 
is not necessarily small. 

In the following sections, we show that by aligning the phase of the coarse and fine solvers, 
it is possible to design parareal algorithms that use multiscale coarse integrators. 

3 Multiscale parareal 

In this section, we introduce the main contribution of this paper - accurate and convergent 
parareal algorithms that use multiscale methods as coarse integrators. Two parareal schemes 
are presented. The first focuses on approximating only the slow variable, while the second 
achieves sup-norm convergence in the state variable, u € M. d . Both methods are based on a 
phase alignment strategy, which can be applied if, for fixed slow variables, the phase is ergodic 
with respect to a circle. Accordingly, we assume that the slow coordinate £ = (£i,■ ■ • ,£d-i) 
is a vector of d — 1 functionally independent slow variables. 

3.1 Multiscale coarse integrator 

For the remainder of this paper, we shall assume that the coarse propagator is a multiscale 
method that only approximates the slow variables. In order to emphasize this point, the 
multiscale coarse integrator will be denoted Alt in place of Cf . Similar to assumption (jl.fi)) . 
we shall assume that 


\SM t x - SM t y\ < (1 + tCHMx) - £(y)|. (3.1) 

The parareal proof of convergence as given in Section 11.21 hinges on the stability assumption 
m, which does not directly involve the exact solution. As a result, as long as m holds, 
the parareal iterations will converge, although not necessarily to the exact solution. However, 
with a multiscale coarse integrator, (13.11) implies that stability only in the slow variables is 
guaranteed. Accordingly, we propose to modify the coarse multiscale integrator by fixing the 
fast variable (a fast phase in the case of HiOsc problems). In terms of slow-fast coordinates, 
the multiscale integrator will be stable in the slow coordinates due to (|3.1I) while stability in 
the fast variable will be enforced by aligning trajectories with respect to a common reference 
phase. In order to achieve this, we assume that one can devise the following local alignment 
algorithm. 

Local alignment: 

Given uo and vq such that £(uo) = £(uo) + A£. 

Let wo = 'I'~ 1 (£( u o)> 4>(vo)) be the point that has the same slow coordinates as uq and the 
same phase as vq. 

Find a point wq such that |u)o — tuol = 0(A£). 
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In other words, the local alignment procedure replaces uq by a new point wq that has the 
same (to order A£) slow coordinates, i.e. £ values, as uq, and approximately the same phase 
as Vo- A trivial solution to the local alignment problem is to set wq := uo- However, this is 
not an adequate strategy that can be used in the next steps of development of our multiscale 
parareal algorithm. 

Notation 3.1. We denote such a local alignment procedure as wq = So(uo;vo)- 

Given a local alignment algorithm Sq, we propose the following modified parareal scheme. 

Algorithm 3.2. 

1. Initialization: (Construct the zero’th iteration approximation) 

Uq =uq and u Q n = n = 1,..., N. 


2. Iterations: k = 1... K 

(a) Parallel fine integrations for n = k ,..., N, 

u F,n = ?HU k ~_ 

(b) Parareal correction: For n = k ,..., N , 

u 0 = u 0 an d U k = So(Aif{Un_i; Up n ) + Up n — So(A4HUn—l'i u F,n)- (3-2) 


In each iteration we first calculate all fine scale integrations. Then, the results of the multi¬ 
scale integrators are aligned with the fine scale ones. In the following, we prove that using 
Algorithm 13.21 all slow variables converge to their limiting value given by the fine scale ap¬ 
proximation. We consider a 1st order multiscale integrator with local phase alignment. 

Theorem 3.3. Let K < N/2 = T/2H. Then, for all k < K, 


sup 

n=0 


£(«n) ~ £ (FuHUq) 


< CH k . 


Proof: We recall the assumption that there exists a diffeomorphism T : u —>• (£(u),<j)(u)) 
such that £ o u(t) are slow while <f> o u(t ) are fast. The variables (£, </>) are only used in the 
analysis but not in the numerical algorithm. 

The main difference with the general analysis described in Section 11.21 is that the bound 
(11.61) is not valid if a multiscale coarse integrator is used. Instead, denoting by 5So(M.tUi] u*) = 
<So(3htti; u*) — So(MtUi',u*), we have 


T o 8S 0 (M t ui\u*) -To 5So{Mtu 2 ; u*) = (<5£, S(f>) , (3.3) 


such that |(5£| < (1+C't)|£(rti)—£( 7 / 2 )| but \5cj)\ = 0(e) which is the accuracy of local alignment. 
Comparing with conventional methods as a coarse integrator and the related estimate (11.121) . 
the slow part is controlled by the local phase alignment in Algorithm 13.21 just like in the non¬ 
singular case, while the rapidly changing phase is incorrect but does not affect the accuracy 
of the slow variables. 
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The slow variables of u k in (13.21) are 


£(“n) = £(So(MHUn-l> u F,n )) + £(. u F,n) ~ £,(<So(MhU 


k-1 


-li u F,n))i 


which is valid with the local alignment. We may thus think of the multiscale integrator 
combined with the local alignment as a coarse integrator with first order accuracy for the 
slow variables. For shorthand, we denote by the combined The 

error in the slow variables is evaluated similarly to (11.81) . 


£(“n) ~ £{FnHU o) 


^(MHUn-l) +€(J 7 HUn_\) - €(MhU*_\) ~ i{^H^{n-l) H Uo) 


k- 


^(MhUu-i) - F(n-\)H u o) + £{SMHF(n-l)H u o) ~ 


+ 


C(5J r H U l ^_ 1 1 ) - (,{5FHF{n-l)HUo) 


Using (13.31) . for every slow variable £, we have that 


H u n) - Ci^nHUo) 


< 


(1 + CH) a<-i) - ^(n-DHUo) + C(E f + H)H £«li) 


Denoting 9 k = (1 + CH) k ~ n (E f + H)~ k H~ k |£(u*) - £ (^nH«o)| 
procedure as in (11.91) . we have for the slow variable, 


and following the same 


sup 

n=0,...,N 


£( u n) - £ (FnHUo) 


< C(Ef + CH) k < CH k . 


3.2 Phase continuity in the coarse and fine scale simulations 

We next consider convergence of the parareal approximation to the exact solutions. The 
main idea is to enforce consistency in the fine scale solutions between neighboring coarse time 
intervals. We may rephrase this problem as the following. 

Forward alignment of step size H : 

Given uq, vo, and u\ = EhUq such that £(«o) — C( v o) = 0(e)- 
Let wq = ^ _1 (£(«o),<Kuo)) and wi = E H w 0 . 

Find a point w\ such that £(u>i) = £(uq) + 0(e) and <f(w i) = 4>i w i) + 0{H 2 ). 


In the problem of forward alignment, if wq is a point with the same slow variable as uq and 
phase as vq, then a forward alignment procedure constructs an order H 2 approximation of 
w\ = J~hWq, the right end point of a coarse interval. See Figure [2]4 for a schematic sketch. 

Notation 3.4. We denote such a forward alignment procedure as w\ = Sjj(ui;uo,vo), where 
A are precomputed parameters to be used in the alignment. 

The forward alignment procedure can be trivially accomplished simply by setting w\ to 
v\EhVq or w\. However, this would require the additional computation of v\ from uq, and so 
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this trivial ’’fix” has a computational cost of sequentially solving the entire system with the 
fine integrator. In practice, for the purpose of parallel in time computations, one needs to 
do so with a computational cost that is lower than running the fine scale solver sequentially. 
Hence, we need to estimate the solution to the given ODE with the given initial condition 
vq by certain simple operations performed on the fine scale solutions already computed in 
parallel. In the following section, we shall describe a forward alignment algorithm for the 
special case of HiOsc ODEs, in which, for fixed slow variables, the fast phase is periodic. The 
method applies only a local exploration by means of minimal additional fine scale computation 
of the solution around uq and u\. In particular, its efficiency is independent of e. 

To summarize, we present the complete multiscale-parareal algorithm. Recall that Mh 
is a multiscale method that only approximates the slow variables. For the fast phase, using 
local and forward alignments, we propagate the needed phase adjustments sequentially along 
with the parareal correction. 


Algorithm 3.5. Full multiscale-parareal algorithm. 

1. Initialization: Construct the zero’th iteration approximation: 

«o = uq and = M.hUu- 1 i n=l,...,N. 

2. Iterations: k = 1...K 

(a) Parallel fine integrations for n = k ,..., N: 


u 


k-1 


r« : = 1 


(b) Header: For n = 0,..., k — 1, set u k = u F ~^. 

(c) Parareal step: Set the initial reference point u* = u k z\ and for n = k,... , N, 

i. Locally align the previous u k 2 Z\ with the current reference point u*: 

u k n Z\ = S 0 (u k n Z\-,u*). 

ii. Align forward to the end of the coarse segment: 

iii. Corrector: 

u k = SoiMnu^] u k F ~') + u k F ~' - S 0 (M H u^z\-,u k F J : ). 

iv. Update the reference point u* = ttjj and repeat. 


Local and forward alignment steps (Step 2(c)i and ii, respectively) create a point ujL 1 at the 
end of each coarse segment according to which all points in the current corrector iteration can 
be aligned. Since the error in each forward alignment is of order H 2 , we find that the overall 
phase is continuous up to a global 0(H ) error. We conclude that following phase alignments, 
the aligned coarse multiscale method provides an globally 0(H ) approximation of both slow 
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and fast variables, i.e., it approximates the solution in the sup norm. 

Example 16.21 demonstrates the effectiveness of the method in a more complicated expand¬ 
ing spiral with a slowly changing frequency. Before describing numerical methods for local 
and forward alignments of HiOsc ODEs, we address the convergence of the algorithm. 


3.3 Convergence of Algorithm 13.51 

Convergence of Algorithm 13.51 in the state variable is obtained in two steps. First, following 
Section m and Theorem 13.31 all slow variables converge to their values obtained by the fine 
scale integrators, 

sup |£(«£) - £(u(nH))\ < C(H k + E f ), 

n=0,...,N 

where C is a constant that is independent of e. In particular, if Ef = O(e), then, following 
0(log(e)) iterations, the error in the slow variables is of order e. As a result, after a few (typ¬ 
ically one or two) iterations, the assumptions underlying forward alignment, that the error in 
the slow variables is of order e holds (more precisely, €(uo) — £(vq) = 0(e)). We may thus think 
of the adopted multiscale combined with the local/forward alignment algorithm as a coarse 
integrator with first order accuracy for all state variables. Hence, the conventional parareal 
proof of convergence described in Section 11.21 holds. More precisely, suppose that, following 
the phase alignment, the set {wq, ... ,wn} is computed and updated in every iteration. The 
following estimates hold for j = 0,..., N, 

\Z(w j )-£(u(jH))\=0(E m ), 

| <f>(w j )-<j>(u(jH))\=0(H), 

\w j -u{jH)\=0(H) + 0{E m ), 

where E m is the error of the aligned multiscale method in approximating the slow variables. 

Assume further the stability properties for the fine and aligned-multiscale coarse propa¬ 
gators m- Therefore, after one parareal iteration, 

sup \u l n - u(nH)\ < C(e~ l E m + Ef). 

n=0,...,N 


After k iterations, and assuming a first order multiscale coarse integrator, E m = O(H), 


sup 

n=0,...,N 


i{u k n ) - t(u(nH)) < C(H k + E f ), 


and 


sup 

n=0,...,N 


u k -u{nH) <C(e 


-i H k 


+ E f ). 


The accuracy of slow variables is improved to by a factor of H per parallel iteration (compare 
with the diverging factor of {e~ l H) k in (11.12(1 ). 

Remark 3.6. In m, Legoll et al propose a multiscale parareal algorithm for stiff ODEs in 
which the fast dynamics is dissipative, i.e., trajectories quickly converge to a lower dimen¬ 
sional manifolds. Unlike the HiOsc case, a naive application of the parareal methodology to 
stiff dissipative systems converges. However, it is not very efficient and suffers from similar 
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difficulties as discussed earlier. To circumvent these difficulties, Legoll et al m suggest a cor¬ 
rection step that allows a consistent approximation of the fast-slow dynamics with parareal. 
This work assumes that the system is split into slow and fast variables, or alternatively, that a 
change of variables that splits the system into such coordinates is given explicitly. Essentially, 
the idea of m is to set the fast component of the multiscale solver with that obtained from 
the fine one. This step may be viewed as a simple alignment method. Indeed, if the (£, (ff) 
coordinates are know, then the same approach can also be applied to the HiOsc case. In 
contrast, the method presented in the following section is seamless in the sense that it does 
not require knowing the slow nor the fast variables. 



Figure 2: Local and forward alignments. (A) At t = 0, given two points uq and vo, we wish 
to approximate wq = 5o(uo; vo) - a point that has the same slow variables as uq and the same 
phase as vo- At t = H, we approximate the point w\ = Fhw o (B) At t = 0, a small 0(e) of 
wq yields at t = H a larger 0(H 2 + e) error. The center of the circles represent w o and w\. 


4 Phase alignment strategies 

In this section, we describe a numerical method for both local and forward alignments as 
defined in the previous section for the special case of HiOsc ODEs in which, for fixed slow 
variables £, the dynamics of the fast phase <f> is periodic. 

In the Algorithm 13.51 vo and uq will correspond to u ^ and u k n ~ 1 , the solutions computed 
at the current and the previous iterations, respectively. The assumption is that vq is the 
more accurate approximation of the solution at the time t = nH , particularly in the phase 
variable. The goal is that from the available information, uq, vo, and u\ := J-huq, we estimate 
v\ := Thv o at t = [n + l)H in order to make correction in the phase of u\. We also emphasize 
that in the subsequent time steps, FhUq is always available because of the prior parallel fine 
integrations. Now wo, as defined in Section [3721 is a point on the same slow coordinates as 
uo but has the same phase as vq. Consequently w\ := is a good estimate of v\. In 

this section, we propose a strategy that move uq to wo, and u\ to a state that is within 
0(e) to wi, without computing Fhw o or FhVq. Our goal is to describe a method that finds 
a point wo such that |rco — Lt>o| = O(e) (local alignment) and a second point w\ such that 
|uq — f&i| = 0(H 2 + e) (forward alignment). 

In addition to Assumption 12.21 we assume the following, 

Assumption 4.1. The fast variable <j> € M and go(£, 4>) is 1-periodic in <j>. 


14 














For fixed £, the time derivative of (j) may depend on the slow variables, i.e., the periodicity 
in time of ffo(£, 0(f)) is of order e and depends on £. Accordingly, it is denoted er(£), where r 
is a smooth, slow function. Note that this does not mean that the oscillation in the original 
state variables are linear because the transformation T is in general nonlinear. 


4.1 At time t move uq closer to w$ 

Assuming that Z( u o) — Z( v o) = 0(e), we may use vq instead of wq, which is not known. Denote 

J(t-,uo,vo) = \F t u 0 - v 0 \ 2 ■ 


We look for the local minima of J(t;uo,vo) closest to t = 0 (by the periodicity assumption, 
such local minima exist) 


0 = J'(t-,u 0 ,v 0 )(t) = 2 (F t u 0 - v 0 ) ■ ^ (F t u 0 - v 0 ) 


2 (Ftu o — vq) ■ 


9'I' 


-i 


dZ 


< + 



To leading order in e, we have 


(T tU0 -v 0 )-(d^~ 1 /d^) = O(e). ( 4 . 1 ) 

In other words, the phase of is close to that of vq, (j){FtUo) = if>o + 0(e) and therefore 
also to the phase of wq. We denote the “first” two local minima 

—eTo + 0(e 2 ) < tg < 0 < < eTo + 0(e 2 ), 

where = r(£o)- Consider 

“'o = F t ±u 0 , 

and the convex combination using these points, 

w 0 = Su 0 = X+Wq + X-Wq, 

with weights A = (A+, A_) independent of e. Thus, equation (14.11) implies that, for any linear 
combination such that A + + A_ = 1, |rco — uo| = 0(e), i.e., wq defined above is a valid choice 
in the local alignment procedure. We define 

Local alignment: 

So{u 0 ]Vo) = X+F.+uq + X-F.-uq. 

c o z 0 


In the numerical implementation of the local alignment So, we use a simple algorithm which 
solves the I 2 minimization of J(t) with a small 0(e) step size and adaptively increasing 
the search domain until an appropriate minimum is found. Denoting the time where the 
minimum is attained by t*, we improve the accuracy by a quadratic interpolation through 
the neighboring points of t*. This algorithm achieves 0(e) accuracy with an efficiency that is 
independent of e. See Appendix [A] for details and [51 j for further references on other efficient 
minimization techniques. 
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4.2 At time t + H move u\ closer to w\ 

We would like to do the same at t = H, i.e., move u\ = FhUq to w\ = FhWq. The main 
difficulty is that we cannot expect that the solution has oscillations of constant periodicities. 
We denote ri = t(£i) = tq + (£1 — £o)(<9r/c0) + 0(H 2 ). In analogy to the procedure at t = 0, 
we find the “first” two minimizers of 


J(t]ui,wi) = \F t ui - wi\ , 


such that 

—eri + 0(e 2 ) < tf < 0 < tf < er\ + 0 (e 2 ). 


Let 


w\ = + A -W 1 , with wf = T t ±u\ = w\ + 0 (e). 


Then, for any constants A + , A_ we have that \w\ — w\\ = 0(e). The problem is that we do not 
know w\ and therefore cannot find tf. One option is to use t± instead and choose weights A± 
that minimize the error. This requires us to relate the t^s and ti’s. 

Denote 

*«o = (Co, 0o) , = d'T'ffUo = (£i, fii ), 

*uo = (%,0o), 

= (Co,0o), Vwi = 'SiFhwq = (Cl, 0i) • 


Without loss of generality, we assume that 0o > 0o and |0o — 0o| < 1. Similarly, assume 
01 > (f\ and |0i — fii | < 1. Then, to leading order in e, 


= (V’o - 0o)eho, = 

to = -(! - 00 + 00)eT 0 , tf = -(1 - 0 i + 0 i)en. 


Next, denote the solution of (£,0) with initial condition (Co,0o) as C(i;Co,0o) and 0(t;Co,0o), 
i.e., (C(0 Co, 0o), 0(0 Co, 0o)) = ^t«o- Using the averaging principle ( 12 . 21 ) . we can write 

£(L Co, 0 o) = i(t) + e 70 /e, C) + 0{e 2 ), 


where £(t) is a slow function that does not depend on the phase and 7 ( 5 , C) is independent of 
e and is r(£)-periodic in s with zero average, Jq^ 7 (s,£)ds = 0 . 

r H 

0(#;Co,0o) = 00 + / [e 1 5'1 (€(*; Co, 0o)) + 32 (C( 0 Co, 0 o))] dt 
Jo 

= 00 +/ [e _1 5i(C(i)) +52(C0))] <& + / si (£(t)h{t/e, €(t))dt + 0(e) 

io JO 

= 0o + U(Co; e) + 0(e), 


for some function F that depends only on £0 and e, but not on the initial phase 0o- In 
particular, we note that 


0Oi) -0o = 0(J0'tuo) -0o = 0(tf; Co, 0o) - 0o = +(Co;e) + C>(e). 


Similarly, 


01 - 00 = 0(TffU O ) - 00 = 0(14; Co, 0o) - 0o = U(Co;e) + 0(e). 
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Hence, 4>(w i) — tp 0 = (f\ — cj>Q = 0(e). In other words, starting at tco instead of uq introduces 
a phase shift that is practically constant. We have then 


tf = Oi - = Oo - 4>o) e 


dr . . 

T 0 + "H7 ' vsl — £o) 


+ 0(eH 2 + e 2 


dr 


= (V>o - </>o)eT 0 + (V’o - </>o)e^ • (6 - Co) + 0{eH 2 + e 2 ) 
= t(j~ + + 0(eH~ + e 2 ), 


where A = -n—% • (Cl - Co) = 0(1). Similarly, 


Ht 0 d£ 


— to T A + 0(eH~ + e 2 ). 


Consider 


Ftf u i - ^-Htftf U 1 + + e 2 ) - + 0(eH~ + e 2 )- 

Expanding around wf 

F t ±u i = uif - j^A + 0(i7 2 + e 2 ), 

for some vector 5 € W l independent of e. Therefore, taking a linear combination A+ + A_ = 1 
and denoting A = (A+,A_), 

uo, vo) = \+F t +ui + \-F t -ui ( 4 . 2 ) 

= (A+u)^ + A—U^ ) H—-5 (A-|_fg" + A—tp ) A + 0(H~ + e 2 ) 

= ici H ——S (A + tQ + A_tg ) A + 0(H~ + e). 

Finally we see that with the choice 


A+ 


~tn 


_ f— 
°0 L 0 


A_ 


t 


+ 

o 



the first order term cancels. Thus, we obtain a second order accurate forward alignment to 
w\. See Figure [3] for the error of (|4.2I) for the simple example of the expanding spiral (11.111) . 

Algorithm 13.51 applies the convex combination (14.21) in forward alignment. Indeed, con¬ 
vergence in the state variable heavily relies on this step because the new point u ^T 1 after the 
forward alignment is assigned as the reference for local alignment in the next coarse interval. 
See Step 2(c)iii. We emphasize that taking S^ as (14.2jl may shift the slow coordinates of the 
resulting u\ from what was computed by the multiscale coarse integrator and assumed to be 
accurate. In the next subsection, we propose a more elaborate method to further improve the 
overall accuracy of the forward alignment step. 


4.3 Improving accuracy in forward alignment 

Here, the idea is that we identify the convex combination with the point which divides the 
trajectory of (11.101) originating from J- f + u i and ending close to F t - u\ by a proportion of A_ 


17 













Figure 3: The error in correcting the phase at the end of one coarse segment for the expanding 
spiral (11.111) . (A) In local alignment, the phase at the end of a coarse segment no is aligned 
with no with an O(e) error. (B) Blue: forward correction with s + , red: backward correction 
with S- and black: a linear combination of shifts using the forward alignment algorithm 
defined in Section |4j With the proposed linear combination, the error is of order H 2 . 


to A + . Since there are two orientations of defined by forward and backward in time 

integrations, the modified convex combination will provide us with two points depending on 
the orientations, and we will choose the one closer to S^u\. 

First, we propose to find the first two local minimizers of 

J(t] Ju+ni, F t -u\) = \J r t(J 7 t +ui ) — T t -u\\ 2 , 

z o z a z o z 0 

such that —eri < < 0 < T+ < er\. Denoting 


t+ + = 4 + A_ r ; and =4 + X-T~, 
we again find the first local minimizers of 

J(t] J-.-ui, J-.++ui) = |J+(J 7 ,-ui) — J r ,++ui| 2 , 

z o z o z o z o 

J{t] Ff-ui, F t +-u\) = u\) — T)+-ui| 2 , 

such that —eri < < 0 and 0 < T+ < eri, and denote them by 

*o“ =t 0 +r;, ^o + = *o + r t- 


(4.3) 


(4.4) 


(4.5) 


With local minimizers of f)4.4|> . the phases between J~ t ++ and and between J~ t +- and 

J 7 t -+ are the same. Now, we define the new weights using t^s in (14.31) and (14.51) by 


A ++ — 


*o \ _ t 0 

— , *+- — 7+3 


f++ 


+++ _ + 
L 0 ''0 


f -1— _ f —I" ’ 

L o L o 


A_= 


+++ _ + 

L o L o 


—) A_ + — -p 


/T-_ f — 

L o L o 


The convex combination (14.21) is now modified as 
^(-Uis-Ucnuo) = A++J r t ++ui + A— F t —u\, S^ 2 (ui-,u 0 ,v 0 ) = X+-J 7 t +-ui + \- + F t -+ui. 
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Here, we note that 

^(S^(ui;uo,v 0 )) = £(iti) + 0(e), f(S^ 2 (ur,u 0 ,v 0 )) = £(ui) +0(e). 

In words, the modified convex combinations (u \; uq. vq) and (tq; uq, Vo) guarantee the 
accuracy of order e in the slow variables of u\. 

Now, we propose to implement the forward alignment Sjj(ui;uo,vo) as follows. 

Forward alignment: 

1. Set the reference point using ()4.2f) . u := \ + F t +u\ + \-F t -u\. 

2. Compute two modified convex combinations with opposite orientations, 

\ ++ F.++u\ + A—Aj— F.+-u\ + A— uJ-.-+u\. 
r 0 l 0 l 0 z 0 

3. Denote by S^(u±] uq, vq) the combination closer to u. 


Remark 4.2. An unperturbed system of the HiOsc system () 1.10 1) . if exists, preserves the slow 
variables but changes the fast variables. Indeed, by denoting F® the fine integrator for the 
unperturbed system, £(J-)°uo) = £(uo) for all t > 0 but (j)(Ffuo) / 4>(uo) for some t > 0. If 
the unperturbed system of ( 11.1011 is explicitly known, one can achieve more accurate local and 
forward alignments by using F° in the minimization of J(t). Unfortunately, the unperturbed 
system is not explicitly known for general systems. Comparison of the parareal solutions using 
F with F° will be presented in Sections m and 16.61 

5 A multiscale integrator based on Poincare-map 

Even though the goal of the multiscale system is a consistent description of only the slow 
variables, in practice, obtaining an explicit expression for the slow variables is often difficult 
or impossible, in particular for high-dimensional systems (see [|9] for an example). In (2], a 
new type of multiscale methods using a Poincare map technique was introduced. This method 
only assumes the existence of slow variables but does not use its explicit form. A novel on- 
the-fly filtering technique achieves high order accuracy. Recall the general two-scale ODE 
(11.101) with initial condition uq E D C M d : 

u = e~ l fi{u) + f 0 (u), u(0)=uq. (5.1) 

By ignoring the lower order perturbation part of the vector field, an unperturbed dynamical 
system is defined. The essential part of the Poincare-map technique is to generate a path whose 
projection on the slow subspace has the correct slow dynamics. To this end, the scheme solves 
both the perturbed and the unperturbed systems from the same initial conditions for short 
time intervals, and compares the resulting trajectories. 

The method relies on the following assumptions regarding the HiOsc dynamics 

Assumption 5.1. The dynamics of the unperturbed equation 

v = e~ 1 f 1 (v), v(0)=v o . (5.2) 

is ergodic with respect to an invariant manifold M(vq). 
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We denote the solution of (15.21) by v(t\v q). 


Assumption 5.2. The invariant manifolds A4(z) is defined by the intersection of the level 
sets of slow variables £i,£ 2 , ••• >£fcj k < d. More precisely, we may identify the invariant 
manifold of v by level sets of the slow variables for u, M(z) = n^ =1 {C € M. d : £j(z) = £j(£)}. 

Hence, the solution u(t) defines a foliation of invariant manifolds M(t) := Ai(u(t)). Note 
that our method only assumes the existence of such £’s but does not require obtaining them. 

Suppose we solve the full equation (15.lj) and the associated unperturbed version (15.21) with 
the same initial condition. Then, it is possible to extract the flow of M(t) from comparison 
of u(t) and v(t) without explicitly knowing the slow variables. The central idea is to locally 
create a path 7 in states space that is transversal to the fast flow. This cut will be defined 
and approximated by a procedure that realizes a Poincare return map along it. We shall look 
for a slow 7 (t), i.e., require that |q| < C such that for any slow variable £, £(7 (t)) = £(u(f)). 
In other words, the effective slow path 7 (t) goes through the same foliation of slow manifolds 
as the exact solution, M.{p/ (f))) = Ai(u(t)). The time derivatives of such effective paths can 
be obtained by extracting the influence of lower order perturbations in the given oscillatory 
equation. Approximating the derivative will require solving the HiOsc system for reduced time 
segments of order e. Since 7 is slow, it can be approximated using macroscopic integrators 
with step size H independent of e. As a result, the overall computational complexity of the 
resulting algorithm is sublinear in e -1 . 

To be consistent with previous notation, we denote by Ft the fine scale approximated 
propagator for the full equation ( 11 . 101 ) . and by Ff the fine scale approximated propagator 
for the unperturbed equation (15.21) in which the low order perturbation is turned off. In 
particular, note that under the dynamics of (15.21) all slow variables are constants of motion. 
Let e < rj < H . A basic Forward-Euler step, depicted in Figure 0]A can be written as 


u n -\-i — FrjUn H — (^Fi^u n F v u n ) . 


(5.3) 


The values of the effective path 7 (t) at t = nH is then identified with u n , n = 0,1, • • • ,N. High 
order approximations of 7 (t) may be obtained by combining several steps and using high-order 
extrapolation. The name ’’Poincare technique” alludes to the fact that 7 (t) is transversal to 
the solution curves of the unperturbed equation. Thus, the full solution induces a Poincare 
return map, which is used to approximate 7 (t). See [ 2 ] for details. 

The bottleneck in the efficiency of the new algorithm is a consequence of small-amplitude, 
high-frequency oscillations in £(u(f)). The accuracy can be improved by sampling the deriva¬ 
tives of a locally smooth average of £, £ instead of the weakly oscillating £. Since we assume 
no explicit knowledge about the slow variables, £ must be computed intrinsically. In m , a 
filtering technique is proposed for the simple case in which the invariant manifold of the un¬ 
perturbed equation is diffeomorphic to a circle, i.e., the unperturbed dynamics is periodic. 
More precisely, we propose to replace (ll.lOp by the filtered equation 


u = e 1 fi(u) + K(t/rj)fo(u,t,e 1 t), t* < t < t* + 7 , (5.4) 


where the filter K r] is C q ([0, r/]) and satisfies the moment condition of the form 



= 0 , 1 , 2 , • • • ,p. 
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The accuracy of (|5.4I) was demonstrated and analyzed in [2]. 




Figure 4: The Poincare map-type technique approximated all slow variables but does not 
require knowing their explicit formulas. (A) A forward-Euler type construction. (B) The 
symmetric Poincare method is first order, but symmetric with respect to the Poincare return 
points. 

5.1 Symmetric Poincare methods 

The simple Forward-Euler step ()5.3|) can be applied in simple situations in which the frequency 
of the fast oscillation is not a slow variable itself (i.e., g± in (12.11) is not a function of £) [2j. 
This restriction can be lifted by generating interpolation points 7 jj* symmetrically described 
as follows. Our idea it to generate and choose interpolation points 7* so that in the state 
space, (277) —1 (7^ — 7AJ approximates implicitly the derivative of £(t) but results in small 
derivative of 4>(t), more precisely, of order rj 2 . The method, originally proposed and analyzed 
in @3], can be described as 

u n+ i = 7* 1 + ^ (7* - 7-i) , (5-5) 

where 

7-1 = rfun, 7* = J^n^Un. 

Convergence of the method is proved in Appendix B. See also [43] for different types of 
Poincare methods. The method ([5. 5 1) defines a propagator, denoted M.h, 

The efficiency of the combined parareal with multiscale Poincare method can be evaluated by 
counting the number of fine steps in each parareal iteration. An effective time segment r is 
denoted the computational cost of a single parareal iteration equivalent to fine integration of 
length r: 

t = H + h fine T] (- — -+ 1 \ (l + ^) (5.6) 

where hfi ne , hp 0 i ncare , and h p h aS e are the microscopic step sizes used for fine, Poincare and 
phase alignment methods, respectively. Overall, with K iterations, the computational speed¬ 
up (assuming maximal parallelization) compared to direct numerical simulation is T/Kt. 
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6 Numerical examples 

6.1 Expanding spiral I 

Consider the following nonlinear equation in the complex plane 

u = ie -1, u| , u| + u\u\~ l , (6-1) 

with the initial value u(0) = 1. The associated unperturbed equation is known to be 

v = ie~ 1 v\v\, u(0) = vo- (6-2) 

As in EH, the dynamics of u(t) can be analyzed by the corresponding system of slow and 
fast variables: 

n = i, m = h 

m = o. ( ‘ } 

We see immediately from Definition 12.11 that £ is a slow variable. The difficulty in the phase 
alignment lies in the singular term in the RHS of <f>. Note that (16.31) is never used in our 
algorithm as £ in (16.3p is only used to show the convergence in the slow variables. 

In this example, we use Algorithm 13.51 ODE45 as a fine integrator and the Poincare 
2nd order multiscale method as a coarse integrator (the Midpoint rule macro-solver and 
ODE45 micro-solver with z-shape construction of 7) to compute the solution. In each micro¬ 
simulation, we solve the filtered equation 

u n (t) = ie _1 u n (t)|u n (i)| + K v (t - i n )w n (i)|« n (t)| _1 , t n < t < t n + 77, (6.4) 

with a C 4 kernel with p = 1 supported on [0, 77] is used. The parameters are specified in 
Table [2] where r]p 0 incare and hp 0 i ncare are parameters for (16.41) . The absolute errors in the 
slow variables and in the state variables as a function of parareal iterations for each different 
e = 10~ 3 and 10“ 4 are shown in Figure [5j In addition, we compare the errors when both local 
and forward alignments are established using the full system ()6.1I) and the unperturbed system 

(16.21) . The parareal solutions using (16.11) in phase alignments have an error of 0(e) in the state 
variables. On the other hand, the phase alignment using (16.2j) shows an error in the state 
variables comparable to that in the slow variables. Indeed, knowing the unperturbed equation 

(16.21) is an advantage. In applications, obtaining an explicit expression for the unperturbed 
equation may not be possible, particularly for nonlinear systems. 


Table 2: Parareal parameters in Example 16.11 


T 

H 

hfine 

V Poincare 

^Poincare 

RelTol, AbsTol (ODE45 parameters) 

2 

1/10 

e/200 

7e 

e/10 

10~ 13 , 10 _ii 
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Figure 5: Expanding spiral I, Example 16.11 (Left) Plot of ||£(-) — £o n fc (-)||i,°°([o,T]) as a 
function of iteration and (Right) the absolute error in the state variables as a function of 
iteration. 


6.2 Expanding spiral with slowly varying fast oscillations 

Consider the following HiOsc example 

x = —27re _1 [l + (1 — az\)z2}y + bx 
y = 27re -1 [l + (1 - azi)z 2 ]x + by 

(6-5) 

Zl = l 

Z2 = -az 2 

where a,b > 0 are constants. Initial conditions are (x,y,zi,Z2)(0) = (1,0,0,1). The solution 
of (16.511 is x(t) = e bt cos [27re -1 (l + e~ at )t \, y(t) = e bt sin [27re _1 (l + e~ at )t ], zi(t) = t, and 
z 2 {t) = e~ at . 

Hence / = x 2 + y 2 , z\ and Z 2 are three slow variables while (x, y) is a linear oscillator with 
expanding amplitude y/l and a slowly changing period e/(l + e~ at ). The example falls under 
the general category of HiOsc systems in which the dynamics of the fast phase slowly evolves 
according to the slow variables. The main difference between this example and Example 16.11 
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is that the derivative of slow variables is not a constant. As a result, the local 0(H 2 ) error 
introduced by the 1st order coarse multiscale integrator is realized. 

The system (16.51) is integrated using the full multiscale Poincare-parareal method, applying 
the corrected phase shift described in Section [4] to ensure convergence in the state variable. 
We stress that the numerical approximation is obtained without using our knowledge that 
the system can be decomposed into the three slow variables /, z\ and z 2 and a fast phase-like 
variable (f> = arctan(y/x). This decomposition and the exact solution are only used in order to 
explain the fast-slow structure in the dynamics and for demonstrating the rate of convergence 
of different variables. 

Figure [TK shows the error in the slow variables as a function of iteration. After a single 
iteration the error in the slow variables drops below e, which is the theoretical limit possible 
with multiscale methods on their own. Figure[7|3 shows the absolute error in the state variable 
of the entire trajectory. Initially, the absolute error is large. This is because the inaccurate 
slow variables create a jump in the phase between coarse time segments. However, after two 
iterations, the phase shift becomes accurate and the method converges to the exact solution. 
Parameters are detailed in Table [4j 


Table 3: Parareal parameters in Example 16.21 


e 

T 

H 

hfine 

VPoincare 

^Poincare 

RelTol, AbsTol(ODE45 parameters) 

10" 3 

2 

0.1 

e/200 

7e 

e/10 

10" i3 , 10 _ii 



Figure 6: Expanding spiral with a = 0.2, b = 0.1, Example 16.21 (A) The error in the slow 
variables, maxj|£j(-) — o u fc (‘)llL°°([o,T])) as a function of iteration. (B) The absolute error 

in the state variables as a function of iteration. 


6.3 Stellar orbits in a galaxy 

The following system is taken from the theory of stellar orbits in a galaxy [391 HO] : 

{ f\ + a 2 ri = er|, 

T 2 + b 2 r 2 = 2erir 2 . 
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Following a change of variables x = [x±,vi,X 2 , V 2] 7 = [ r i> if r i/ a , r 2 , 'jjf^/b]' 1 ' and after a 
rescaling of time, t = er, the system can be written in the following form 

x\ = e~ l av i 

vi =-e~ l axi+x 2 2 /a 

X 2 = e V 2 

v 2 = -e~ l X 2 + 2xi x 2 /b. 

Initial conditions are (xi,vi,X2,v 2 )(0) = (1,0,1, 0). Resonance of oscillatory modes take effect 
in the lower order term when a = 2 and 6=1. Using the algorithm proposed in [5], three 
independent slow variables are identified as 

= xj + vj, £ 2 = x 2 2 + vl , £ 3 = x 1 xl + 2vix 2 v 2 -xivl. (6.7) 

The example falls under the category of HiOsc systems in which two stiff harmonic oscillators 
are coupled. The local 0(H 2 ) error introduced by the 1st order coarse multiscale integrator 
is realized. 

The system (| 6 . 6 D is integrated using Algorithm 13.51 applying the global alignment algo¬ 
rithm described in Section |4] to ensure convergence in the state variable. The Poincare method 
is used as a coarse solver, using the trajectory of (16.61) as the flow T and the one of (16.61) 
without a lower order perturbation as the flow J 70 . We stress that the numerical approxima¬ 
tion is obtained without using our knowledge that the system can be decomposed into the 
three slow variables £ 1 , £2 and £3 and a fast phase-like variable <fi. This decomposition are 
only used in order to explain the fast-slow structure in the dynamics. 

Figure [7] shows the absolute error in the state variable of the entire trajectory. Initially, 
the absolute error is large because the inaccurate slow variables create a jump in the phase 
between coarse time segments. However, after four iterations the error in the state variable 
drops below e, which is the theoretical limit possible with multiscale methods on their own. 
Parameters are detailed in Table U The fine integrator is ODE45 method, and the coarse 
integrator is the Poincare 2nd order multiscale method. A C 3 kernel with p = 1 is used for 
the filtered equation. 


Table 4: Parareal parameters in Example 16.31 


e 

T 

H 

hfine 

Poincare 

^Poincare 

RelTol, AbsTol(ODE45 parameters) 

10” 4 

14 

0.5 

e /100 

20e 

e /10 

10~ ia , 10" 11 


6.4 Non-linear oscillators 

Consider the following example of a Voltera-Lotka oscillator with slowly varying frequency 
and amplitude 

x = e _1 x(l — zy) 

V = e~ l zy{x - 1) 
z = 0.2x 
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Figure 7: Stellar orbits in a galaxy, Example 16.31 The absolute error in the state variables as 
a function of iteration is depicted. 


Initial conditions are (x,y,z)( 0) = (1,2.9,1). For fixed z, ( x,y) is a Voltera-Lotka oscillator 
whose period is of order e. The period and amplitude of ( x , y ) depend on a parameter z, 
which is given by the time integral of x. As a result, z is a slow variable. It is easily verified 
that the first integral of the oscillator is also slow, 

I = X - log(x) + y - log (y)/z 

Again, we stress that the slow variables are only used in order to demonstrate the results of 
the method. They are not used in the numerical approximation. In addition, Figure [8]A shows 
the level set of the slow variable, {u € M 3 : I(u) = I(x(t n ), y(t n ), z(t n ))} , projected onto x-y 
plane. In contrast to the previous examples, the level set of the slow variable I is not a circle. 
As a result, J{t) may have several local minima and we need to find the first local minima 
which is close to the global minimum of J within a few periods. Parameters are given in Table 
[H The fine integrator is ODE45 method, and the coarse integrator is the Poincare 2nd order 
multiscale method. A C 3 kernel with p = 1 is used for the filtered equation. 

Table 5: Parareal parameters in Example 16.41 


e 

T 

H 

hfine 

VPoincare 

^Poincare 

RelTol, AbsTol(ODE45 parameters) 

10” 3 

10 

1/2 

e/200 

30e 

e/10 

10~ 13 , 10 -1U 


6.5 Passage through resonance 

One of the fundamental assumptions underling multiscale approaches such as Poincare and 
other methods is a spectral gap in the spectrum of the Jacobian of the equations of motion. 
When this assumption fails, for example due to a temporary passage through resonance, the 
assumption 13.11 may not hold close the resonance and typical multiscale methods fail. How¬ 
ever, the applicability of the multiscale algorithms can be extended by the parareal approach 
described above, i.e., by resolving all scales of the dynamics - both the slow and the fast. 
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Figure 8: A Lotka-Voltera oscillator with slowly varying frequency and amplitude, Exam¬ 
ple [631 (A) The level sets of the slow variable I projected onto x-y plane. (B) The absolute 
error in the state variables as a function of iteration. 


We consider the following example. 

x = —2ne~ 1 f(z)y + 0.5 sin(z)x 
y = 2ne~ 1 f(z)x 
z = 1 


where 

f(z) = tanh (50(2: — 4.5)). 

Initial conditions are (1,0,0). In words, f(z) changes smoothly from -1 to 1, vanishing at 
2 ; = 4.5. Hence, the frequency of oscillation undergoes fast oscillations with varying frequency, 
except close to t = 4.5. At this time, f(z) vanishes and the system is no longer highly 
oscillatory. More precisely, trajectories go through a transition layer. Its width in this example 
is of order e. The two slow variables are I = x 2 + y 2 and 2 . 

Figure [9j\ shows the values of the state variables with e = 10~ 4 . Due to the resonance, 
the Poincare method fails to capture the correct evolution of the slow variables when crossing 
the singular point t = 4.5. However, combining with parareal, the fine solution of parareal 
integrates the equation across the resonance and allows the multiscale method to proceed 
beyond the singularity. In Figure [9)3, the absolute error in the state variable does not decrease 
with iterations because the accuracy of phase alignment relies on the scale separation which 
does not exists near t = 4.5. We show, however, that the convergence in the state variable can 
be achieved with a slight modification of Algorithm 13.51 Figure [9]0 is obtained by skipping 
phase alignments, Step 2(c)i and ii, and replacing Step (c)iii with the naive correction (11.211 
in the interval near t = 4.5. 

Remark 6.1. We note that this example goes beyond the scope for which convergence is proven 
in Section [H Indeed, the purpose of the example is to demonstrate that the applicability of 
the multiscale-parareal coupling using alignments may be wider than proven here. 
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Table 6: Parareal parameters in Example 16.51 


e 

T 

H 

hfine 

?iPoincare 

^Poincare 

RelTol, AbsTol(ODE45 parameters) 

1(T 4 

7 

1/4 

e/200 

15e 

e/10 

10~ 13 , 10 _ii 



t k k 


Figure 9: Passage through resonance, Example 16.51 (A) The solution of x(t) and y(t) with 
t € [4,5] and e = 10~ 4 . The frequency function vanishes at t = 4.5 and solutions lose their 
highly oscillatory nature. (B) The absolute errors of both state (circles) and slow (crosses) 
variables as a function of iterations with phase alignment at all the time. (C) The absolute 
errors with phase alignment turned off near t = 4.5. Convergence in the state variable is 
achieved. 


6.6 The Fermi-Pasta-Ulam (FPU) problem 

We consider a chain of 2k springs on a line, connected with alternating soft nonlinear and stiff 
linear springs with both ends fixed. This problem has been used as a benchmark for testing 
the long-time performance of geometric integrators m- See also mini so] and references 
therein for related recent work. The model is derived from the following Hamiltonian: 

^ 2 k ^ k k 

H(p, q ) = + 4 e ~ 2 X^ 2i ~ 92*-i) 2 + 5^(©i+i - Q2i) 4 - 

i=l i —1 i —1 

Using the change of variables given in j5], the equations of motion for the system can be 
written as 


Vi = Ui , 
xs = 

. _ _ ( 6 . 8 ) 

Ui = -{Vi- £Xi - yi -1 - exi-iY + (y i+ i - ex i+ i - yi - exi) A , 

Vi = ~e~ l Xi + (yi - exi - yj_i - exj__i) 3 + (y i+1 - ex i+1 - yi - exi) 3 . 

Since the ends are fixed, yo = xq = ]Jk+i = a^fc+i = 0. 

In this section, we endeavor solving (16.81) in the 0(e~ 1 ) time scale using the proposed 
methods. The aim is to expose the limit of the various different algorithmic components of 
the proposed methods. Clearly, the time scale in which we compute the solutions of (16.81) is 
out of the scope of the analysis that we presented earlier. 
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System (|6.8I) is solved in M 8 with k = 2, and it admits seven slow variables — they are 
the total energies of the stiff springs, fi = xf + vf for i = 1,2, the relative phases between 
the stiff springs, <j> = X\X 2 + V\V 2 , and all the degrees of freedom which are related to the soft 
springs: yi and tq, i = 1,2. The nontrivial energy transfer and the relative phase take place 
in the very long e _1 time scale. 

In Figure [TOlA. we present the maximum errors in the state variable in a long time interval 
[0, T] = [0,e _1 /2], with e = 10~ 3 . The system (|6.8D is integrated using Algorithm 13.51 
applying the global alignment algorithm described in Section 0] to achieve convergence in the 
state variable. The results are computed with the initial conditions y\ = x\ = j /2 = ^2 = lj 
(u\, v\i U 2 , V 2 ) = (0,1.2,1,0) and with the parameters given in Tabled The Verlet method 
was used as the fine integrator and the Poincare 2nd order method (Verlet macro-solver and 
ODE45 micro-solver) as the multiscale integrator. We further point out that the errors will 
decrease as e becomes smaller. 

In Figure flOtl. we present the maximum errors computed using the same set of conditions 
and parameters, except that the phase alignments are computed by the numerical solutions of 
the modified equation: 

Vi = 0 • Ui, 

Xi = e~ l Vi, 

< (6.9) 

Ui = 0 ■ {-(yi - exi - yi -1 - exi- 1) 6 + (y i+ 1 - ex i+1 - yi - ex*) 3 } > 

Vi = -e~ l Xi + 0 ■ {(yi-exi- y t -\ - ex*_ 1) 3 + (y i+1 - ex i+ i - yi - ex*) 3 } , 

with yo = xq = yk+i = Xk+i = 0. This is a slightly more general type of “unperturbed” 
equations as those defined in Section 5. This approach may be thought of as phase alignments 
with some slow variables constrained. A detailed discussion about such topic is out of the 
scope of this paper, and shall be addressed in a future work. 

Remark 6.2. We observe that in this study, phase alignments by solutions of (16.911 yield 
significantly superior results. Indeed, our experience indicated that for the problems that can 
be solved by the Poincare method described in Section [H it is generally better to use the 
so-called “unperturbed” equations for phase alignments, as the slow variables are constrained 
and would not deteriorate the performance of the phase alignment steps. However, we chose 
to present a more general setup in this paper, avoiding any specific choice of coarse/multiscale 
integrators in describing the algorithms. 


Table 7: Parareal parameters in Example 16.61 


e 

T 

H 

hfine 

V Poincare 

^Poincare 

RelTol, AbsTol(ODE45 parameters) 

nr 3 

e-72 

1/4 

e /20 

10e 

e /20 

10~7 10” iU 
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Figure 10: The Fermi-Pasta-Ulam problem, Example 16.61 The absolute error in the state 
variables as a function of iteration is depicted. (A) Phase alignments with the full system 
(16.81) . (B) Phase alignments with the unperturbed equation explained in the paragraph. 


7 Summary 

The paper describes two approaches to incorporate multiscale integrators as coarse integrators 
in parareal methods. The first, presented in Section 18. 11 approximates all the slow variables. 
However, the numerical approximation of the state variables does not converge to the true 
solution u(nH). This parareal-multiscale combination has several advantages compared to 
other multiscale schemes. 

• It offers increased stability and is less sensitive to the choice of parameters. Intuitively, 
the parareal iterations can ’’fix” errors incurred by the inexact multiscale scheme. 

• It offers increased accuracy. In fact, the accuracy of slow variables may be smaller than 
0(e), which is a theoretical limit for Poincare and other multiscale methods that are 
based on averaging or homogenization principles. 

• It may be applied for systems with moderate scale separation. Most multiscale methods 
are more efficient than conventional, non-multiscale schemes if the separation in scale is 
large enough, i.e., if e is sufficiently small. However, they typically become less efficient 
or unstable at intermediate values of e. 

• The method may be used in situations in which the dynamics looses its multiscale 
structure in a short transition layer, for example, due to passage through resonance, see 
Example 16.51 

The second approach, presented in Section 13.21 computes convergent approximation to 
all state variables in the system. This algorithm requires the phase alignment procedure, 
described in Section |4] in addition to the steps needed in our first algorithm. We prove that 
the accuracy of the scheme in the sup norm after K iterations is of order e~ l H K + Ef, where 
H is the a coarse step size. In particular, the number of iterations to achieve a given error 
tolerance is logarithmic in e -1 . 

The computational cost of the method can be divided into two contributions. The first 
is the cost of the fine integrator invoked at each parareal iteration. With K iterations its 
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contribution to the overall cost is proportional to KHe~ l . The second contribution comes from 
the overhead of coarse multiscale integrators and phase alignment. While this contribution is 
independent of e, it grows linearly with the number of coarse step sizes, H -1 . Hence, there is a 
trade off in choosing H. With a large scale separation e< 1, the first contribution dominates 
and, assuming maximal parallelization is available, it is advantageous to use a relatively small 
H, even if the multiscale method allows larger steps. The two contributions balance if one 
takes H = y/e, which implies a computational cost of order Ke ~ 1 / 2 . In contrast, parallel 
methods using conventional integrators will require at least 0(e ~ x ) steps. 

The multiscale-parareal coupling strategy using alignments is not limited to the proposed 
particular implementation involving a single-frequency fast variable. In the multi-frequency 
case, different alignment algorithms may be possible. For example, if the slow-fast system 
is given explicitly, one could simply set the fast variable in the coarse integrator to coincide 
with the fine one. 
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A Appendix: Adaptive search algorithm 

We explain our implementation for the local alignment 5o- This algorithm adaptively searches 
local minima of the functional J(t) by adjusting the step size and the computational inter¬ 
val. The first part of the algorithm numerically solves the I 2 minimization problem of J(t) 
using quadratic interpolation. The second part computes a convex combination for the local 
alignment So (uq ; vq ). 

For simplicity, we assume that the fine solver applies a numerical method with step size 
h = 0(e). Also, let r] = r] p h aS e > 0 denote a parameter that is larger than at least one period. 

It is assumed that 77 is of order e, but the size of 77 is not explicitly known. Our algorithm 
adaptively finds the size of 77 and locally integrates the HiOsc ODE in short time segments of 
length [—77,77], Thus, the computational cost of each alignment procedure is independent of e. 

Recall that at t = 0, given two points uq and vq such that |£(rto) — £(uo)| = 0(e), local 
alignment approximates a new point w 0 such that £(u>o) = £(«o) and cj>(wo) = In 

addition, first forward and backward alignment times t± are sought, in which \T f ± uq — uo| are 
local minima. 

In the following, denote by [x\ the integer value of x. 

Algorithm A.l. Local alignment Sq(uo',vq) 

1. Adaptive search for first two local minima of J(f) associated with positive(-f-) and 
negative(-) orientations. 

(a) Set 77 = e, h = e/100, Tol = 00, u prev = 00, and I± = {}. 

(b) While Tol > e/100 or I± + {} 

i. Forward/Backward local integration: Compute uo t i = TihUo, i € {0,..., [±77 /h\}. 
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ii. Compute the local minimum of J(i) = |uo,« — v o\ 2 , * € {0,, [±r]/h]}. 

Denote it I± = (ii±, i 2 ±, • • • , i r ±) ■ 

iii. If I± ^ set Tol — \up rev ^o,q±l’ 'Uprev — ,u±? nnd h — h/ 2. 

Else, set = 2 p and h = e/100. 

End while. 

(c) Quadratic interpolation: For each of the first two indices i\±, improve the min¬ 
imization of J using a quadratic interpolation. Let p±{t) denote the polynomial 
such that 

P±(h± + j) = u 0 ,i 1± +j, j = -1,0,1. 

Denote by t ^ the minimum of | p±(t) — uq|. 


2. Let 


and define 


A 4- = 


t 0 +t 0 


A_ = 


to 


t 0 + t 0 


<So(uo',vo) = X+^+uq + A-J^-uq. 


(A.l) 


In (IA.1I) . we implement the fine integrator T with step size h determined in step 1(b). If t±/h 
is not integer, then T is to be evaluated using quadratic interpolation. 


B Appendix: Convergence of the symmetric Poincare method 

To prove convergence of the symmetric Poincare method algorithm described in Section 15.11 
we use a diffeomorphism T : u —» (£(u), (f>(u)) given in (|2.1D . 

ff = 5o(£,0), £(0)=£(«o), / B ^ 

\0 = e _ V(£) + S 2 (£, 0 ), 0(0 )=ct>(u 0 ), 

where go(^,4>) and 52(^1 0) are 1-periodic in 0. We stress that the variables (£,0) are only 
used in the analysis but not in the numerical algorithm. In Section [2l the dynamics of the 
slow variables can be approximated by an averaged equation of (12.21) . 


I = HO, F(0 = J <*,(£, 0)c% 

m = £(uq). 

Recall that the construction given in Section T5.II 

TJ 

* 1 ( * * \ 

Un +1 = 7-1 + 2 ^ (7i - 7 - 1 ) • 

To simplify the calculation, we generate 7 * using the symmetric shape (z-shape) which is 
centered at 7^. See Figure The formulas for 7* 1; 7 q and 7* are thus of the form 

7-1 = FpF-pUn, 7q = U n , 7 * = 

We then prove the following theorem. 
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Theorem B.l. The Poincare force estimator defined by 


V — _(F 0 T — T°T \ 

' ~ 2 ? y v r -v' rT i • r v' r ~ r i) 

satisfies the following estimates in the (£, 0) coordinate: 

f{Vu n ) = J t i{t n ) + E f , \4>(Vu n )\ =o(j£ 


(B.3) 


where Ep denotes the error of the filtered equation in approximating an averaged equation. 


Proof. Consider the unperturbed system associated with (IB.111 . 

f£ = 0, £(o) = £o, 

10= 751 (0) 0(°) = 0o- 


(B.4) 


When T is corresponding to the filtered equation (15.41) . it is proved in [2| that F± v u n is essen¬ 
tially close to averaged £(t n ±r]), respectively. We denote by (£ + , 0 + ) the slow-fast coordinate 
of the point 7 ^ i.e., (£ + ,0 + ) = (£( 7 *), 0(7*))- Similarly, (£ _ ,0 _ ) = (£( 7 * 1 ), 0(7-i))- Then 
£ + and 0 + satisfy 

U + =fo + J^F{fi F {t))dfi 

\0 + = 0o + Jo ^ 9 i{£ F {t))dt + /J 7 G(f F (t))dt - f 5 i(£ + ), 

where f F (t) corresponds to the solution of (IB.21) forward in time and G is the averaged 52 (C) 0) 
due to the filtered equation. On the other hand, for 7 * x , 

fr=&-J7«P, 

\(j)- = 00 - /o' ±gi(f B {t))dt - f Q ri G(f B {t))dt + 2gi(C), 

where f B (t) corresponds to the backward in time solution of ()B.2I) . 

For the slow variables, evaluating the force by yields 

£ + -r = [" F(t F (t))dt+ [' F{t B {t))dt 

Jo Jo 

which approximates the evolution of £ over the interval [— 7 , 7 ]. For the fast variable, 

0 + -0“ = [ -9i(€ F (t))dt - -5i(£ + ) + [ -9i(Z B (t))dt - -51 (£ _ ) 

Jo e e Jo e e 

' -v-' 

h 


+r G(£ F (t))dt- r G(£ B (t))dt. 

Jo Jo 


I 2 


By considering £ F (t) = £0 + Jo F(fi F )ds and f B (t) = £0 — jJ F(fi R )ds, we can compute the 
Taylor series to estimate el\, developed at £q. 
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eh = ^ gi (^0 + F (£ F (s))ds^ dt-rig! (^,0 +F(£, F (t))d?J^ 

+ | J 9i (to ~ F(^ B (s))ds\ dt - ngi ^ F(£ B (t))dt^J j 

= y* yy 9 i(t;o) T F(z F (s))ds^ 

~{Jo yj^^ofFi^is^ydt-vJ^Vg^ofF^imtyOir, 3 ). 

Here, the terms gg\ (£q) canceled. Due to the symmetric structure, expanding both ( s ) and 

£ b (s) about s = 0 using Taylor series, all 0(g 2 ) terms vanish and we have 

eh = dis¬ 
similarly, one can show that 

h = 0( V 2 ). 

Therefore, there exist nonnegative constants Ci and C 2 such that 

3 

|0(7i) - 0 ( 7 * 1 )| = I0 + - 0“l < Ci^- + C 2 V 2 . (B.5) 

This completes the proof. □ 

It is important to use an appropriate sequence of interpolating points in the state space, 7 ^, 
as it directly impacts on the accuracy of Poincare method. With a parameter r/ = Ce, (IB.3|) 
shows that the force estimator V generates two interpolating points to estimate the evolution 
of the slow variables over 2g with an 0(e) disagreement in the fast variable. This results in 
stable and accurate approximations. However, we remark that the force estimation using the 
algorithm (15.31) sometimes introduces an 0(1) difference in the fast variable and thus shifts 
the slow variables. 
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